////////////////////////////////////////////////////////////////////////
//MCMC_EdgeOp.cpp
//
// edge operator: Log
//				  Canny
//
// Author: Zhuowen Tu
////////////////////////////////////////////////////////////////////////
// color image

#include "mcEdgeOp.h"
#include "opencv2/highgui/highgui.hpp"

#define NUM_DIRECTIONS 8

////////////////////////////////////////////////////////////////////////
Raster<double> NonMaximalSuppress(const Raster<double> &Mag,
						const Raster<double> &Ix,
						const Raster<double> &Iy,
						const double alfa,
						const bool	bSpecial)
{
    int		x,y;
    double	g;               /* center pixel gradient magnitude value */
    double	dx, dy;
    double	g1, g2;          /* interpolated gradient magnitude values */
    double	ga, gb, gc, gd;  /* nearest pixels' gradient magnitude values */
    double	ux, uy;          /* components of edge normal vector */

	Raster<double> I;

	// initialize output image
	I.SetDimension(Mag.rows(),Mag.cols());
	I.Zero();

	double mag_min, mag_max;
	mag_min = Mag.GetMin();
	mag_max = Mag.GetMax();
	// magnitude threshold
	double threshold=alfa*(mag_max-mag_min)+mag_min;

    /* For each pixel (except those on the image boundary), check if
       the gradient magnitude is a local (3x3 neighbourhood) maxima in
       the gradient orientation: */
    for (x=2; x<Mag.rows()-2; x++)
	{
		for (y=2; y<Mag.cols()-2; y++)
		{
			dx = Ix.Data(x,y);
			dy = Iy.Data(x,y);

			if (bSpecial)
			{
				// not possibly being an edge point
				if (fabs(dx) <= EBSHLON && fabs(dy) <= EBSHLON)
					continue;

				g = Mag.Data(x,y);
				if ( ((g>Mag.Data(x-1,y) && g>Mag.Data(x+1,y)) ||
					 (g>Mag.Data(x,y-1) && g>Mag.Data(x,y+1)) ) &&
					g > threshold)
					I(x,y) = g;
				else
					I(x,y) = 0;
			}
			else
			{
				if (fabs(dy)>fabs(dx))
				{
					ux = fabs (dx)/fabs (dy);
					uy = 1.0;

					gb = Mag.Data(x,y-1);
					gd = Mag.Data(x,y+1);

					if (Sign(dx)==Sign(dy))
					{
						ga = Mag.Data(x-1,y-1);
						gc = Mag.Data(x+1,y+1);
					}
					else
					{
						ga = Mag.Data(x+1,y-1);
						gc = Mag.Data(x-1,y+1);
					}
					g1 = (ux*gb)+(uy-ux)*ga;
					g2 = (ux*gd)+(uy-ux)*gc;
				}
				else
				{

					uy = fabs(dy)/fabs(dx);
					ux = 1.0;

					gb = Mag.Data(x+1,y);
					gd = Mag.Data(x-1,y);

					if (Sign(dx)==Sign(dy))
					{
						ga = Mag.Data(x+1,y+1);
						gc = Mag.Data(x-1,y-1);
					}
					else
					{
						ga = Mag.Data(x+1,y-1);
						gc = Mag.Data(x-1,y+1);
					}
					g1 = (uy*gb)+(ux-uy)*ga;
					g2 = (uy*gd)+(ux-uy)*gc;
				}

				g = Mag.Data(x,y);

				if (g>g1 && g>g2 && g>threshold)
					I(x,y) = g;
				else
					I(x,y) = 0;
			}
		}
    }
	return I;
}


////////////////////////////////////////////////////////////////////////
MCMC_Pos & MCMC_Pos::operator=(const MCMC_Pos &a)
{
	if (&a != this)
	{
		i = a.i;
		j = a.j;
	}

	return *this;
}

////////////////////////////////////////////////////////////////////////
MCMC_Edge& MCMC_Edge::operator=(const MCMC_Edge &a)
{
	if (&a != this)
		*((LinkList<MCMC_Pos> *) this) = *((LinkList<MCMC_Pos> *) (&a));

	return *this;
}

void MCMC_Edge::AddAnEdgePoint(const int j, const int i)
{
	MCMC_Pos pos(i,j);

	Insert(pos);
}


////////////////////////////////////////////////////////////////////////
MCMC_EdgeList::MCMC_EdgeList()
{
	p_edge_list = NULL;
	num_edge = 0;
}

MCMC_EdgeList::~MCMC_EdgeList()
{
	if (p_edge_list != NULL)
		delete []p_edge_list;
	p_edge_list = NULL;
	num_edge = 0;
}

void MCMC_EdgeList::SetEdgeNum(const int num)
{
	if (p_edge_list != NULL)
		delete []p_edge_list;
	p_edge_list = NULL;
	num_edge = 0;

	num_edge = num;
	p_edge_list = new MCMC_Edge[num_edge];
}

MCMC_Edge* MCMC_EdgeList::GetFirst(void)
{
	move = 0;

	if (move < num_edge)
		return &(p_edge_list[move]);
	else
		return NULL;
}

MCMC_Edge* MCMC_EdgeList::GetNext(void)
{
	move++;

	if (move < num_edge)
		return &(p_edge_list[move]);
	else
		return NULL;
}

void MCMC_EdgeList::InsertEdgePoint(const int index, const int j, const int i)
{
	if (p_edge_list!=NULL && index>=0 && index<num_edge)
		p_edge_list[index].AddAnEdgePoint(j,i);
}


///////////////////////////////////////////////////////////////////////
MCMC_Corner::MCMC_Corner()
{
	k = 0;
	i = 0;
	j = 0;
}

MCMC_Corner::~MCMC_Corner()
{
}

void MCMC_Corner::Create(const int num, const int n, const int m, mcmcVector<double> &vc_angle1)
{
	assert(num >= 0);

	k = num;
	j = n;
	i = m;
	vc_angle.SetRowVector(k);

	for (int t=0; t<k; t++)
		vc_angle(t) = vc_angle1(t);
}


MCMC_CornerList::MCMC_CornerList()
{
	num_corner = 0;
	p_corners = NULL;
}

MCMC_CornerList::~MCMC_CornerList()
{
	Free();
}

void MCMC_CornerList::Free(void)
{
	num_corner = 0;
	if (p_corners != NULL)
		delete []p_corners;
	p_corners = NULL;
}

void MCMC_CornerList::Create(const int num)
{
	// free allocated space
	Free();

	assert(num >= 0);
	num_corner = num;
	p_corners = new MCMC_Corner[num_corner];
}

MCMC_Corner* MCMC_CornerList::GetCorner(const int i)
{
	if (i>=0 && i<num_corner)
		return &(p_corners[i]);
	else
		return NULL;
}

MCMC_Corner* MCMC_CornerList::GetFirstCorner(void)
{
	move = 0;
	return GetCorner(move);
}

MCMC_Corner* MCMC_CornerList::GetNextCorner(void)
{
	move++;
	return GetCorner(move);
}

//generate a matrix mask where pixles labeled as '1' are corner points
void MCMC_CornerList::GetnerateCornerMask(void)
{
	// 1: generate matrix mask
	ras_corner_mask.InitValue(0);

	MCMC_Corner *pcorner;

	pcorner = GetFirstCorner();
	while (pcorner!=NULL)
	{
		ras_corner_mask(pcorner->J(),pcorner->I()) = 1;
		pcorner = GetNextCorner();
	}
}

///////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////
MCMC_EdgeOp::MCMC_EdgeOp()
	:Raster<double>()
{
	dx[0]=dy[4]=dy[5]=dy[6]=dx[6]=dx[7]=-1;
    dx[1]=dx[5]=dy[3]=dy[7]=0;
    dx[2]=dx[3]=dx[4]=dy[0]=dy[1]=dy[2]=1;
}

MCMC_EdgeOp & MCMC_EdgeOp::operator=(const Raster<uchar> &I)
{
	CopyRaster(*((Raster<double>*) this), *((Raster<uchar>*) &I));

	return *this;
}

MCMC_EdgeOp & MCMC_EdgeOp::operator=(const Raster<double> &I)
{
	CopyRaster(*((Raster<double>*) this), *((Raster<double>*) &I));

	return *this;
}

// edge points: 255
Raster<double> MCMC_EdgeOp::Canny(double sigma, double alfa, Raster<double> *pIx,Raster<double> *pIy, bool ifspecial)
{
	// generate kernel along x direction
	mcmcMatrix<double>	h1,h2;
	Raster<double>	Ix,Iy;
	Raster<double>	I;

	h1.SetGaussian(1,sigma,0);
	h2.SetGaussian(1,sigma,1);

	// x direction
	Ix = ConvHor(h2);
	Ix = Ix.ConvVer(h1.T());
	// y direction
	Iy = ConvHor(h1);
	Iy = Iy.ConvVer(h2.T());

	// get the magnitude of the first dirivitive of image
	Raster<double>	Mag;
	int			i,j;

	Mag.SetDimension(rows(),cols());

	for (j=0; j<rows(); j++)
		for (i=0; i<cols(); i++)
			Mag(j,i) = sqrt(pow(Ix(j,i),2.0)+pow(Iy(j,i),2.0));

	// nomaximal suppression
	I = NonMaximalSuppress(Mag,Ix,Iy,alfa,ifspecial);

	// copy gradient
	if (pIx != NULL)
		*pIx = Ix;
	if (pIy != NULL)
		*pIy = Iy;

	return I;
}




/////////////////////////////////////////////
// color image edge operators

// construction & destruction
MCMC_EdgeOp_Color::MCMC_EdgeOp_Color()
	:Raster_Color<double>()
{
	dx[0]=dy[4]=dy[5]=dy[6]=dx[6]=dx[7]=-1;
    dx[1]=dx[5]=dy[3]=dy[7]=0;
    dx[2]=dx[3]=dx[4]=dy[0]=dy[1]=dy[2]=1;
}

/////////////////////////////////////////////
// operators
 MCMC_EdgeOp_Color & MCMC_EdgeOp_Color::operator=(const MCMC_EdgeOp_Color &I)
{
	*((Raster_Color<double>*) this) = *((Raster_Color<double>*) &I);
	return *this;
}


MCMC_EdgeOp_Color& MCMC_EdgeOp_Color::operator=(const Raster_Color<uchar> &I)
{
	CopyRasterColor(*((Raster_Color<double>*) this), *((Raster_Color<uchar>*) &I));
	return *this;
}

MCMC_EdgeOp_Color& MCMC_EdgeOp_Color::operator=(const Raster_Color<double> &I)
{
	CopyRasterColor(*((Raster_Color<double>*) this), *((Raster_Color<double>*) &I));
	return *this;
}

/////////////////////////////////////////////
// operations

double MCMC_EdgeOp_Color::calcPercentThreshold(cv::Mat mImage, float dPercent, int nBins) {
	if (mImage.channels() != 1) return -1;

	int nPixels = cols()*rows();
	int channels[] = { 0 };
	int histSize[] = { nBins };
	double min, max;
	minMaxIdx(mImage, &min, &max);

	float range[] = { min, max };
	float sumN = 0, dMarker = nPixels*dPercent;
	const float* ranges[] = { range };
	cv::Mat hist;

	calcHist(&mImage, 1, channels, cv::Mat(), hist, 1, histSize, ranges);
	for (int i = 0; i < nBins; i++) {
		sumN += hist.at<float>(i);
		if (sumN >= dMarker) {
			return ((max - min) / nBins * i);
		}
	}
	return -1;
}

Raster<double> MCMC_EdgeOp_Color::Canny(double sigma,double alfa,
								  Raster<double> *pIx,Raster<double> *pIy,
								  const bool	bSpecial)
{
	if (!TUs_CANNY) { // opencv version of Canny
		cv::Mat mRGB, mGray, mBlurred, mEdges, Ixy, Ix, Iy;
		float threshHigh, threshLow, median, alpha = 0.33;
		Raster<double> I;

		mRGB = RGB;
		mRGB.convertTo(mRGB, CV_8UC3);
		cv::cvtColor(mRGB, mGray, CV_RGB2GRAY);
		cv::GaussianBlur(mGray, mBlurred, cv::Size(0, 0), sigma);

		//	Calculate the gradient along two directions
		Ixy.create(rows(), cols(), CV_64F);
		Sobel(mBlurred, Ix, CV_64F, 1, 0, 3);
		Sobel(mBlurred, Iy, CV_64F, 0, 1, 3);
		for (int y = 0; y < rows(); y++)
		for (int x = 0; x < cols(); x++)
			Ixy.at<double>(y, x) = sqrt(pow(Ix.at<double>(y, x), 2) + pow(Iy.at<double>(y, x), 2));

		//	Canny edge detection

		cv::Mat mI[3], mE[3];
		cv::split(mRGB, mI);
		for (int i = 0; i < 3; i++) {
			/*
			cv::Scalar mu, sigma;
			cv::meanStdDev(mI[i], mu, sigma);
			threshLow = mu.val[0] - sigma.val[0];
			threshHigh = mu.val[0] + sigma.val[0];
			*/

			// Otsu's threshold: when image consists of distinct foreground and background
			cv::Mat mTest;
			threshHigh = cv::threshold(mI[i], mTest, 0, 255, CV_THRESH_BINARY + CV_THRESH_OTSU);
			threshLow = threshHigh *0.5;

			/*
			median = calcPercentThreshold(Ixy, 0.5, 150);
			threshLow = 0.66 * median;
			threshHigh = 2 * 0.66 * median;
			*/
			cv::Canny(mI[i], mE[i], threshLow, threshHigh, 3, true);
		}
		mEdges = cv::Mat(mRGB.rows, mRGB.cols, CV_8U);
		for (int i = 0; i < mRGB.rows; i++)
		for (int j = 0; j < mRGB.cols; j++) {
			mEdges.at<uchar>(i, j) = max(max(mE[0].at<uchar>(i, j), mE[1].at<uchar>(i, j)), mE[2].at<uchar>(i, j));
		}

		//cv::Canny(mBlurred, mEdges, threshLow, threshHigh, 3, true);

		I.SetDimension(rows(), cols());
		I.data.setTo(cv::Scalar::all(0));
		for (int i = 0; i < rows(); i++)
		for (int j = 0; j < cols(); j++)
		if (mEdges.at<uchar>(i, j) > 5)
			I.data.at<double>(i, j) = Ixy.at<double>(i, j);


		Raster<double> temp;
		temp.SetDimension(rows(), cols());
		temp.data = Ix.clone();
		*pIx = temp;
		temp.data = Iy.clone();
		*pIy = temp;
		return I;
	}
	else { // Tu's version of Canny
		// generate kernel along x direction
		mcmcMatrix<double>	h1, h2;
		Raster<double>	IxR, IyR, IxG, IyG, IxB, IyB, Ix, Iy;
		Raster<double>	I;

		h1.SetGaussian(1, sigma, 0);
		h2.SetGaussian(1, sigma, 1);

		CreateLUVspace();

		// Jiadi Yang
		cv::Mat mRGB[3];
		cv::split(RGB, mRGB);

		Raster<double> rRGB[3];
		rRGB[0].SetDimension(RGB.rows, RGB.cols);
		rRGB[1].SetDimension(RGB.rows, RGB.cols);
		rRGB[2].SetDimension(RGB.rows, RGB.cols);
		rRGB[0].data = mRGB[0];
		rRGB[1].data = mRGB[1];
		rRGB[2].data = mRGB[2];

		// color in red
		// x direction
		IxR = rRGB[0].ConvHor(h2);
		IxR = IxR.ConvVer(h1.T());
		// y direction
		IyR = rRGB[0].ConvHor(h1);
		IyR = IyR.ConvVer(h2.T());

		// color in green
		// x direction
		IxG = rRGB[1].ConvHor(h2);
		IxG = IxG.ConvVer(h1.T());
		// y direction
		IyG = rRGB[1].ConvHor(h1);
		IyG = IyG.ConvVer(h2.T());

		// color in blue
		// x direction
		IxB = rRGB[2].ConvHor(h2);
		IxB = IxB.ConvVer(h1.T());
		// y direction
		IyB = rRGB[2].ConvHor(h1);
		IyB = IyB.ConvVer(h2.T());

		// Jiadi Yang
		// get the magnitude of the first dirivitive of image
		Raster<double>	Mag;
		int			i, j;
		double			p, q, t, theta;

		Mag.SetDimension(rows(), cols());
		Ix.SetDimension(rows(), cols());
		Iy.SetDimension(rows(), cols());

		for (j = 0; j < rows(); j++)
		for (i = 0; i < cols(); i++)
		{
			p = IxR.Data(j, i)*IxR.Data(j, i) +
				IxG.Data(j, i)*IxG.Data(j, i) +
				IxB.Data(j, i)*IxB.Data(j, i);

			q = IyR.Data(j, i)*IyR.Data(j, i) +
				IyG.Data(j, i)*IyG.Data(j, i) +
				IyB.Data(j, i)*IyB.Data(j, i);

			t = IxR.Data(j, i)*IyR.Data(j, i) +
				IxG.Data(j, i)*IyG.Data(j, i) +
				IxB.Data(j, i)*IyB.Data(j, i);

			double lamda = (p + q) / 2.0 + sqrt((p + q)*(p + q) - 4 * (p*q - t*t)) / 2;
			Mag(j, i) = sqrt(fabs(lamda));
			if (p != q)
			{
				theta = atan2(lamda - p, t);
				Ix(j, i) = Mag(j, i)*cos(theta);
				Iy(j, i) = Mag(j, i)*sin(theta);
			}
			else
			{
				theta = PI / 2;
				Ix(j, i) = Mag(j, i);
				Iy(j, i) = 0;
			}
		}

		// nomaximal suppression
		I = NonMaximalSuppress(Mag, Ix, Iy, alfa, bSpecial);

		// copy gradient
		if (pIx != NULL)
			*pIx = Ix;
		if (pIy != NULL)
			*pIy = Iy;
		return I;
	}
	
}

void MCMC_EdgeOp_Color::LabelJunctions(Raster<int> &ras_edgemap)
{
	int i,j,k,prev;
	int ii,jj,count;
	int w,h;

	w = ras_edgemap.cols();
	h = ras_edgemap.rows();

	for(i=0;i<w;i++)
		for(j=0;j<h;j++)
			if(ras_edgemap(j,i) != 0) ras_edgemap(j,i)=1;

	for(i=1; i<(w-1); i++)
	{
		for(j=1; j<(h-1); j++)
		{
			if(ras_edgemap(j,i))
			{
				prev=(ras_edgemap(j-1,i-1)!=0);
				count=0;
				ii=0;
				jj= -1;
				for(k=0;k<8;k++)
				{
					if(k==1) ii=1;
					if(k==2) jj=0;
					if(k==3) jj=1;
					if(k==4) ii=0;
					if(k==5) ii= -1;
					if(k==6) jj=0;
					if(k==7) jj= -1;
					if((ras_edgemap(j+jj,i+ii) != 0)!=prev) count++;
					prev=(ras_edgemap(j+jj,i+ii)!=0);
				}
				/* actual junction or isolated endpt */
				if((count>4) || (count==2))
					ras_edgemap(j,i) = 2;
			}
		}
	}
}

void MCMC_EdgeOp_Color::ile_propagate(Raster<int> &ras_edgemap, const int x, const int y, const int d)
//int x,y;        /* position to propagate from */
//int d;          /* direction */
{
    int le,ue,dd,i,cl,si,sj;
	int	w,h;

	w = ras_edgemap.cols();
	h = ras_edgemap.rows();

    if(d&1) {   /* odd, so three neighbours */
        le= -1;
        ue=1;
    } else {    /* even, so five neighbours */
        le= -2;
        ue=2;
    }
    for(i=le;i<=ue;i++)
	{
        dd=(d+i+NUM_DIRECTIONS)%NUM_DIRECTIONS;
        si=x+dx[dd];
        sj=y+dy[dd];
		if((si>=0)&&(si<w)&&(sj>=0)&&(sj<h)) {
        	if(ras_edgemap(sj,si) == 2) return;
		}
    } /* not blocked */

    cl = ras_edgemap(y,x);
    for(i=le;i<=ue;i++)
	{
        dd=(d+i+NUM_DIRECTIONS)%NUM_DIRECTIONS;
        si=x+dx[dd];
        sj=y+dy[dd];
		if((si>=0)&&(si<w)&&(sj>=0)&&(sj<h))
		{
     	   	if(ras_edgemap(sj,si) == 1)
			{
       	     	ras_edgemap(sj,si) = cl;
       	     	ile_propagate(ras_edgemap,si,sj,dd);
       	 	}
		}
    }
}


// Label number on edge pixels that form connected line
void MCMC_EdgeOp_Color::LabelNumberOnEdgePixels(Raster<int> &ras_edgemap)
{
	int	w,h,counte,countj,i,j,ii,jj;
	int numsegments,count;
	int si,sj,done,d,pass,found;

	w = ras_edgemap.cols();
	h = ras_edgemap.rows();

	/* do a little bit of error checking */

	for(i=1; i<(w-1); i++)
	{
		for(j=1; j<(h-1); j++)
		{
			if(ras_edgemap(j,i) == 1)
			{	/* edge pixel */
				counte=0;
				countj=0;
				for(ii=-1; ii<=1; ii++)
				{
					for(jj=-1; jj<=1; jj++)
					{
						if((ii==0) && (jj==0)) continue;
						if(ras_edgemap(j+jj,i+ii) == 1) counte++;
						if(ras_edgemap(j+jj,i+ii) == 2) countj++;
					}
				}
			}
		}
	}

	/* junctions are 2, edges are 1 */

	numsegments=0;
	done=0;
	while(!done)
	{
		done=1;
		found=0;
		for(i=1;(i<(w-1))&&(!found);i++)
		{
			for(j=1;(j<(h-1))&&(!found);j++)
			{
				if(ras_edgemap(j,i) == 1)
				{
					/* have to check neighbourhood for junctions */
					found=0;
					count=0;
					pass=1;
					si=i;
					sj=j;
					for(ii=-1; ii<=1; ii++)
					{
						for(jj=-1; jj<=1; jj++)
						{
							if((ii==0) && (jj==0)) continue;
							if(ras_edgemap(j+jj,i+ii) == 2) pass=0; /* dont use */
							if(ras_edgemap(j+jj,i+ii) == 1) count++;
						}
					}
					if((pass)&&(count==2)) found=1;
				}
			}
		}
		if(found)
		{
			ras_edgemap(sj,si) = numsegments+3;
			for(d=0; d<NUM_DIRECTIONS; d++)
				ile_propagate(ras_edgemap,si,sj,d);
			done=0;
			numsegments++;
		}
	}

	/* for closed, no-junctions regions, start anywhere */
    done=0;
    while(!done)
	{
        done=1;
        found=0;
        for(i=1;(i<(w-1))&&(!found);i++)
		{
            for(j=1;(j<(h-1))&&(!found);j++)
			{
				if (ras_edgemap(j,i) == 1)
				{
                    /* have to check neighbourhood for junctions */
                    found=0;
                    count=0;
                    pass=1;
                    si=i;
                    sj=j;
                    for(ii=-1; ii<=1; ii++)
					{
                        for(jj=-1; jj<=1; jj++)
						{
                            if((ii==0) && (jj==0)) continue;
                            if(ras_edgemap(j+jj,i+ii) == 2) pass=0; /* use */
                            if(ras_edgemap(j+jj,i+ii) == 1) count++;
                        }
                    }
                    if((pass==1)&&((count==1)||(count==2))) found=1;
                }
            }
        }
        if(found)
		{	/* texture line or int edge */
			ras_edgemap(sj,si) = numsegments+3;
            for(d=0; d<NUM_DIRECTIONS; d++)
                ile_propagate(ras_edgemap,si,sj,d);
            done=0;
            numsegments++;
        }
    }

	/* check that no 1 pixels remain */
	for(i=1;i<(w-1);i++)
	{
		for(j=1;j<(h-1);j++)
		{
			if(ras_edgemap(j,i) == 1)
			{
/*
 * Modification by Brendan McCane to include isolated pixels as an edge
 * and edges with only two pixels in them.
 */
                ras_edgemap(j,i) = numsegments+3;
                for (ii=-1; ii<=1; ii++)
                for (jj=-1; jj<=1; jj++)
                {
                    if ((ii==0) && (jj==0)) continue;
                    if (ras_edgemap(j+jj,i+ii) == 1)
                        ras_edgemap(j+jj,i+ii) = ras_edgemap(j,i);
                }
                numsegments++;
			}
		}
	}
}


void MCMC_EdgeOp_Color::FormEdgeList(MCMC_EdgeList &edge_list, Raster<int> &ras_numbered_edgemap)
{
	int	w,h,num_edges,i,j;

	w = ras_numbered_edgemap.cols();
	h = ras_numbered_edgemap.rows();

	num_edges=0;
	for(i=0;i<w;i++)
		for(j=0;j<h;j++)
			if(ras_numbered_edgemap(j,i)>num_edges) num_edges=ras_numbered_edgemap(j,i);

	num_edges -= 2;	/* for junctions and edges */

	int	nedge;

	edge_list.SetEdgeNum(num_edges);

	for(i=0; i<w; i++)
	{
    	for(j=0; j<h; j++)
		{
			nedge = ras_numbered_edgemap(j,i)-3;
			edge_list.InsertEdgePoint(nedge,j,i);
        }
	}
}


void MCMC_EdgeOp_Color::ExtractEdges(MCMC_EdgeList	&edge_list, double sigma, Raster<int> *pras_result)
{
	Raster<int>		ras_edgemap;
	Raster<double>	ras;
	int				i,j;


	// apply canny edge operator
	ras = Canny(sigma,0.05);

	// get corresponding edge map
	ras_edgemap.SetDimension(ras.rows(),ras.cols());
	for (j=0; j<ras.rows(); j++)
		for (i=0; i<ras.cols(); i++)
		{
			if (ras(j,i)>0.0)
				ras_edgemap(j,i) = 1;
			else
				ras_edgemap(j,i) = 0;
		}
	// label edge junctions
	LabelJunctions(ras_edgemap);
	// label pixel number
	LabelNumberOnEdgePixels(ras_edgemap);
	// extract edge list
	FormEdgeList(edge_list,ras_edgemap);

	if (pras_result != NULL)
	{
		// for output as demon
		pras_result->SetDimension(ras_edgemap.rows(),ras_edgemap.cols());
		pras_result->InitValue(0);

		MCMC_Edge	*p_one_edge;
		MCMC_Pos	pos;
		bool		bGo;

		p_one_edge = edge_list.GetFirst();
		while (p_one_edge != NULL)
		{
			bGo = p_one_edge->GetFirst(pos);
			while (bGo)
			{
				(*pras_result)(pos.j,pos.i) = 1;

				bGo = p_one_edge->GetNext(pos);
			}

			p_one_edge = edge_list.GetNext();
		}
	}
}

void MCMC_EdgeOp_Color::RidgeEdgeSeeds(const double sigma, mcmcMatrix<uchar> &mx_ridge, mcmcMatrix<uchar> *pmx_edges)
{
	Raster<double>	ras,ras_Gx,ras_Gy;
	mcmcMatrix<uchar>	mx_edges;
	int				i,j,m,n,ii,jj,j_c,i_c;
	double			angle1,angle2;

	// apply canny edge operator
	ras = Canny(sigma,0.05,&ras_Gx,&ras_Gy);

	// get corresponding edge map
	mx_edges.SetDimension(ras.rows(),ras.cols());
	for (j=0; j<rows(); j++)
		for (i=0; i<cols(); i++)
		{
			if (ras(j,i)>10.0)
				mx_edges(j,i) = 1;
			else
				mx_edges(j,i) = 0;
		}

	mx_ridge.SetDimension(rows(), cols());
	mx_ridge.InitValue(0);
	for (j=0; j<rows(); j++)
		for (i=0; i<cols(); i++)
		{
			if (mx_edges(j,i)!=0)
			{
				angle1 = AngleO2PI(atan2(ras_Gy(j,i), ras_Gx(j,i)));
				for (m=-8; m<=8; m++)
					for (n=-8; n<=8; n++)
					if (m*m+n*n<=64)
					{

						jj = j+m;
						ii = i+n;
						if (mx_edges.Valid(jj,ii) && mx_edges(jj,ii)!=0)
						{
							angle2 = AngleO2PI(atan2(ras_Gy(jj,ii), ras_Gx(jj,ii)));
							if (DistanceBtwnOrientation(angle1,angle2) < 20.0)
							{
								j_c = IntClose(((double)(j+jj))/2.0);
								i_c = IntClose(((double)(i+ii))/2.0);
								if (mx_edges(j_c,i_c)==0)
									mx_ridge(j_c,i_c)=1;
							}
						}
					}
			}
		}
	if (pmx_edges != NULL)
		(*pmx_edges) = mx_edges;
}



#define NUMANGLES		12
#define	NUMBD			4
#define VOTE_THRESHOLD	2.5

int MCMC_EdgeOp_Color::DiscoverDirections(mcmcVector<double> &vc_votes, mcmcVector<int> &vc_index)
{
	int		inspected,kept,i;
	int		temp;
	double	temp_vote;

	inspected = vc_votes.size();
	kept = 0;

	while (inspected > 0)
	{
		temp_vote = 0.0;
		temp = -1;
		for (i=0; i<vc_votes.size(); i++)
		{
			if (vc_votes(i)<VOTE_THRESHOLD && vc_votes(i)>=0.0)
			{
				vc_votes(i) = -1.0;
				inspected--;
			}
			else if (vc_votes(i)>temp_vote && vc_votes(i)>=VOTE_THRESHOLD)
			{
				temp_vote = vc_votes(i);
				temp = i;
			}
		}
		if (temp>=0)
		{
			for (i=0; i<kept; i++)
				if (abs(temp-vc_index(i))<=1 || abs(temp-vc_index(i))==NUMANGLES-1)
					break;
			if (i>=kept)
				vc_index(kept++) = temp;
			vc_votes(temp) = -1.0;
			inspected--;
		}
	}
	return kept;
}

// classify detected corners and generate branch angles
void MCMC_EdgeOp_Color::CornerClassification(Raster<double> &ras_Ix, Raster<double> &ras_Iy,
								 Raster<uchar> &ras_corner, Raster<uchar> &ras_edgemap,
								 MCMC_CornerList &corner_list)
{
	mcmcVector<double>	vc_votes(NUMANGLES),vc_angle(NUMANGLES);
	int				i,j,m,n,index,num_corner,k;
	double			angle1,angle2;
	MCMC_Corner		*pcorner;
	mcmcVector<int>		vc_index(NUMANGLES);


	// Discover how many number of corners having at least two branches is.
	num_corner = 0;
	for (j=NUMBD; j<ras_Ix.rows()-NUMBD; j++)
		for (i=NUMBD; i<ras_Ix.cols()-NUMBD; i++)
		{

			// find the corner points
			if (ras_corner(j,i)==1)
			{
				// clean the votes
				vc_votes.InitValue(0.0);

				// find the edge points
				for (n=j-NUMBD; n<=j+NUMBD; n++)
					for (m=i-NUMBD; m<=i+NUMBD; m++)
						if (ras_edgemap(n,m) == 1)
						{
							// gredient angle
							angle1 = atan2(ras_Iy(n,m),ras_Ix(n,m));
							if (angle1 < 0)
								angle1 += PI;
							angle2 = atan2((float)(n-j),(float)(m-i));
							angle1 = angle1*Sign(angle2);
							if (angle1 < 0)
								angle1 += 2*PI;


							index = IntClose(angle1/(2.0*PI/NUMANGLES));
							index = index%NUMANGLES;

							vc_votes(index) = vc_votes(index)+1;
							vc_votes((index+1)%NUMANGLES) = vc_votes((index+1)%NUMANGLES)+0.5;
							vc_votes((index+NUMANGLES-1)%NUMANGLES) = vc_votes((index+NUMANGLES-1)%NUMANGLES)+0.5;
						}

				k = DiscoverDirections(vc_votes, vc_index);
				if (k>=2 && k<=5)
					num_corner++;
			}
		}

	// allocate space
	corner_list.Create(num_corner);
	corner_list.ras_corner_mask.SetDimension(ras_Ix.rows(),ras_Iy.cols());
	corner_list.ras_corner_mask.InitValue(0);

	num_corner = 0;
	for (j=NUMBD; j<ras_Ix.rows()-NUMBD; j++)
		for (i=NUMBD; i<ras_Ix.cols()-NUMBD; i++)
		{

			// find the corner points
			if (ras_corner(j,i)==1)
			{
				// clean the votes
				vc_votes.InitValue(0.0);

				// find the edge points
				for (n=j-NUMBD; n<=j+NUMBD; n++)
					for (m=i-NUMBD; m<=i+NUMBD; m++)
						if (ras_edgemap(n,m) == 1)
						{
							// gredient angle
							angle1 = atan2(ras_Iy(n,m),ras_Ix(n,m));
							if (angle1 < 0)
								angle1 += PI;
							angle2 = atan2((float)(n-j),(float)(m-i));
							angle1 = angle1*Sign(angle2);
							if (angle1 < 0)
								angle1 += 2*PI;

							index = IntClose(angle1/(2.0*PI/NUMANGLES));
							index = index%NUMANGLES;

							vc_votes(index) = vc_votes(index)+1;
							vc_votes((index+1)%NUMANGLES) = vc_votes((index+1)%NUMANGLES)+0.5;
							vc_votes((index+NUMANGLES-1)%NUMANGLES) = vc_votes((index+NUMANGLES-1)%NUMANGLES)+0.5;
						}
				k = DiscoverDirections(vc_votes, vc_index);
				if (k>=2 && k<=5)
				{
					for (n=0; n<k; n++)
						vc_angle(n) = (2.0*PI/NUMANGLES)*vc_index(n);
					pcorner	= corner_list.GetCorner(num_corner);
					if (pcorner != NULL)
						pcorner->Create(k,j,i,vc_angle);

					num_corner++;
				}
			}
		}
	corner_list.GetnerateCornerMask();
}


void MCMC_EdgeOp_Color::CornerDetection(double sigma,Raster<uchar> &ras_corner, Raster<uchar> &ras_edgemap, MCMC_CornerList &corner_list)
{
	Raster<double>	ras_Ix,ras_Iy;
	Raster<int>		ras_output;

	// canny operator
	Canny(sigma,0.2,&ras_Ix,&ras_Iy);
	// classify corners
	CornerClassification(ras_Ix,ras_Iy,ras_corner,ras_edgemap,corner_list);
}


///////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////

MCMC_OpAnEdge & MCMC_OpAnEdge::operator=(MCMC_OpAnEdge &a)
{
	if (this != &a)
	{
		vc_points = a.vc_points;
		list_pixels = a.list_pixels;
	}

	return *this;
}

void MCMC_OpAnEdge::ExpendTo(const int num_new_points)
{
	if (num_new_points <= vc_points.size())
		return;
	
	SimpleMatrix<MCMC_Pos>	vc_temp;

	vc_temp.SetDimension(1, vc_points.size());

	for (int k=0; k<vc_points.size(); k++)
		vc_temp(k) = vc_points(k);

	vc_points.SetDimension(1, num_new_points);
	for (int k=0; k<vc_temp.size(); k++)
		vc_points(k) = vc_temp(k);
}


void MCMC_OpAnEdge::CopyFromPixelList(void)
{
	vc_points.SetDimension(1, list_pixels.Num());

	bool	bgo;
	int		i,j,k;

	k = 0;
	bgo = list_pixels.GetFirst(j, i);
	while (bgo)
	{
		vc_points(k).i = i;
		vc_points(k).j = j;
		k++;
		bgo = list_pixels.GetNext(j, i);
	}
}

void MCMC_OpAnEdge::CopyToPixelList(void)
{
	list_pixels.Empty();

	for (int k=0; k<vc_points.size(); k++)
		list_pixels.Insert(vc_points(k).j, vc_points(k).i);
}

void MCMC_OpAnEdge::CopyFrom(MCMC_OpAnEdge &a, int num)
{
	vc_points.SetDimension(1, num);
	for (int i=0; i<num; i++)
		vc_points(i) = a.vc_points(i);

	CopyToPixelList();
}

void MCMC_OpEdges::ToImage(Raster<int> &ras_edges, bool bone)
{
	int		k,m;

	if (!bone)
		for (k=0; k<NumEdges(); k++)
		{
			assert(Edge(k).vc_points.size() == Edge(k).list_pixels.Num());

			for (m=0; m<Edge(k).NumPoints(); m++)
				if (ras_edges.Valid(Edge(k).j(m), Edge(k).i(m)))
					ras_edges(Edge(k).j(m), Edge(k).i(m)) = k+3;
		}
	else
		for (k=0; k<NumEdges(); k++)
		{
			assert(Edge(k).vc_points.size() == Edge(k).list_pixels.Num());

			for (m=0; m<Edge(k).NumPoints(); m++)
				if (ras_edges.Valid(Edge(k).j(m), Edge(k).i(m)))
					ras_edges(Edge(k).j(m), Edge(k).i(m)) = 1;
		}
}


//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
MCMC_EdgeLinking::MCMC_EdgeLinking()
{
	dx[0]=dy[4]=dy[5]=dy[6]=dx[6]=dx[7]=-1;
    dx[1]=dx[5]=dy[3]=dy[7]=0;
    dx[2]=dx[3]=dx[4]=dy[0]=dy[1]=dy[2]=1;
}

MCMC_EdgeLinking::~MCMC_EdgeLinking()
{
}

void MCMC_EdgeLinking::ObtainAtomicRegions(Raster<double> &ras_canny, Raster<double> &Ix, Raster<double> &Iy, double scale)
{
	// angle
	ras_gradient.SetDimension(ras_canny.rows(), ras_canny.cols());
	ras_angle.SetDimension(ras_gradient.rows(), ras_gradient.cols());
	
	ras_edges.SetDimension(ras_gradient.rows(), ras_gradient.cols());
	
	int i,j;
	
	for (j=0; j<ras_canny.rows(); j++)
		for (i=0; i<ras_canny.cols(); i++)
		{
			ras_angle(j, i) = AngleO2PI(atan2(Ix(j, i), Iy(j, i)));
			ras_gradient(j, i) = min(ras_canny(j, i)/15.0*184, 255.0);

			if (ras_canny(j, i) > 0.0)
				ras_edges(j, i) = 1;
			else
				ras_edges(j, i) = 0;
		}

	Raster<uchar> ras_output;

	CopyRaster(ras_output, ras_gradient);
		
	LinkEdges(0, 25, 1);

	ras_output = ras_edges.ToVisibleRaster();

	// thin edges
	ThinEdges();

	ras_output = ras_edges.ToVisibleRaster();


}

 

/*---------------------- Link up Zero Crossings ---------------------------*/

/* Scan the image and link adjacent non-zero points into lists.
   This is applied only to the first band of the image.
   Note that the image is modified as scanned edge points are zeroed. 
*/

void MCMC_EdgeLinking::LinkEdges(const double th_low, const double th_high, const double min_len1)
{
	threshold_low_linking = th_low;
	threshold_high_linking = th_high;
	min_len = min_len1;

	int i,j,m,n;
	MCMC_PixelSet		tracked_edge_points;
	bool				bgo;
	LinkList<double>	list_gradient;
	double				gradient;

	ras_edges.SetDimension(rows(), cols());
	ras_edges.InitValue(0);

	for (j=0; j<rows(); j++)
		for (i=0; i<cols(); i++)
			if (ras_gradient.Data(j, i) > threshold_high_linking)
			{
				/* Start tracking an edge from the given position.  As each point is added
				   to an edge, zero its pixel in the image.  After the end of the edge
				   is reached, then reverse points in the array and start tracking from
				   the beginning in the other direction.  Finally, the new edge is added to
				   the VEdges structure.
				*/
				tracked_edge_points.Empty();
				list_gradient.Empty();

				// add this starting point
				tracked_edge_points.Insert(j, i);
				list_gradient.Insert(ras_gradient(j, i));
				ras_gradient(j, i) = 0;
				
				// Track to end of connected points in first direction.
				m = j;
				n = i;
				while (NextPoint(&m, &n))
				{
					tracked_edge_points.Insert(m, n);
					list_gradient.Insert(ras_gradient(m, n));
					ras_gradient(m, n) = 0;
				}

			    // Remove points from top of stack that are below high threshold.
				while (tracked_edge_points.GetLast(m, n))
				{
					assert(list_gradient.GetLast(gradient));

					if (gradient < threshold_high_linking)
					{
						tracked_edge_points.DeleteTheLastPixel();
						list_gradient.DeleteTheLast();
					}
					else
						break;
				}
				
				// Start tracking from original point in opposite direction.
				assert(tracked_edge_points.GetLast(m, n));
				while (NextPoint(&m, &n))
				{
					tracked_edge_points.InsertToHead(m, n);
					list_gradient.InsertToHead(ras_gradient(m,n));

					ras_gradient(m, n) = 0;
				}

			    // Remove points from top of stack that are below high threshold.
				while (tracked_edge_points.GetFirst(j, i))
				{
					assert(list_gradient.GetFirst(gradient));

					if (gradient < threshold_high_linking)
					{
						tracked_edge_points.DeleteTheFirstPixel();
						list_gradient.DeleteTheFirst();
					}
					else
						break;
				}
				
				bgo = tracked_edge_points.GetFirst(m, n);
				while (bgo)
				{
					ras_edges(m, n) = 1;
					bgo = tracked_edge_points.GetNext(m, n);
				}
			}
}

static int roff[8] = {1, 0, -1, 0, 1, 1, -1, -1};
static int coff[8] = {0, 1, 0, -1, 1, -1, -1, 1};
/* Given the location of an image point, find the best edge continuation
     from this point to one of its neighbors.  Return the position
     of the new point. */
#define GAP 1
bool MCMC_EdgeLinking::NextPoint(int *rp, int *cp)
{
    int r, c, next_max=0, nextr, nextc;
	int ii, jj;

    /* These arrays give the row and column offsets for each of the eight
       neighbors of a point, starting with those that are 4-connected. */

    /* Check each of the neighbors for a pixel above the low threshold. */
	for (ii=-GAP; ii<=GAP; ii++)
	{
		for (jj=-GAP; jj<=GAP; jj++)
		if ((ii==0)&&(jj==0)) continue;
	    else
		{
			r = *rp + ii;
			c = *cp + jj;
			if (r >= 0  &&  c >= 0  &&  r < rows()  &&  c < cols())
				if ((ras_gradient(r, c)  >=  threshold_low_linking)&&
					(!next_max))
				{
					nextr = ii;
					nextc = jj;
					next_max = IntClose(ras_gradient(r, c));
				}
		}
    }
	if (next_max==0)
		return false;
	else 
	{
		*rp += nextr;
		*cp += nextc;
		return true;
	}
}


/* junction_filter is of form 012 			*/
/*						      3 4 			*/
/*							  576 			*/
/* And file we use file junction_filter2 	*/
/* This has the affect of skeletonisation	*/

void MCMC_EdgeLinking::ThinEdges(void)
{
	int i,j,done;
	int w,h;
	int numchanges;
	int n;
	int minx,maxx,miny,maxy;
	int jmap[256]={1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, 
				   1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 0, 0, 
				   1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 
				   1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 0, 1, 0, 1, 
				   1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
				   1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
				   1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 
				   0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 
				   1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 
				   1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 
				   0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 
				   0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1, 
				   1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 
				   1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 
				   0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1, 
				   0, 0, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1 };

	w = cols();
	h = rows();

	/* gives crude edges - pass through the junction filter */
	done=0;	
	numchanges=0;
	minx=1;
	maxx=w-2;
	miny=1;
	maxy=h-2;
	while(!done) {
		done=1;
		for(i=minx;i<=maxx;i++) {
       		for(j=miny;j<=maxy;j++) {
				if(ras_edges(j,i)) {
					n=(ras_edges(j-1,i-1)!=0);
					n+=2*(ras_edges(j-1,i)!=0);
					n+=4*(ras_edges(j-1,i+1)!=0);
					n+=8*(ras_edges(j,i+1)!=0);
					n+=16*(ras_edges(j+1,i+1)!=0);
					n+=32*(ras_edges(j+1,i)!=0);
					n+=64*(ras_edges(j+1,i-1)!=0);
					n+=128*(ras_edges(j,i-1)!=0);
					if(!jmap[n]) {
						ras_edges(j,i)=0;
						done=0;
						numchanges++;
					}
				}
			}
		}
	}
}


/* These functions go through a nominally thin-edged image and */
/* set all the junctions to a value of 2                       */

void MCMC_EdgeLinking::LabelJunctions(void)
{
	int i,j,k,prev;
	int ii,jj,count;
	int w,h;

	w = cols();
	h = rows();

	
	// to start the non-zero point in ras_edges much be 1

	for(i=1;i<(w-1);i++) {
		for(j=1;j<(h-1);j++) {
			if(ras_edges(j,i)) {
				prev=(ras_edges(j-1,i-1)!=0);
				count=0;
				ii=0;
				jj= -1;
				for(k=0;k<8;k++) {
					if(k==1) ii=1;
					if(k==2) jj=0;
					if(k==3) jj=1;
					if(k==4) ii=0;
					if(k==5) ii= -1;
					if(k==6) jj=0;
					if(k==7) jj= -1;
					if((ras_edges(j+jj,i+ii)!=0)!=prev) count++;
					prev=(ras_edges(j+jj,i+ii)!=0);
				}
				/* actual junction or isolated endpt */
				if((count>4)||(count==2)) {
					ras_edges(j,i)=2;
				}
			}
		}
	}
}

void MCMC_EdgeLinking::IlePropagate(Raster<int> &ras_edgemap, int x, int y, int d) 
//im;		/* map of labels */
// w,h;
//int x,y;        /* position to propagate from */
//int d;          /* direction */
{
    int le,ue,dd,i,cl,si,sj;
	int	w,h;

	w = ras_edgemap.cols();
	h = ras_edgemap.rows();

    if(d&1) {   /* odd, so three neighbours */
        le= -1;
        ue=1;
    } else {    /* even, so five neighbours */
        le= -2;
        ue=2;
    }
    for(i=le;i<=ue;i++)
	{
        dd=(d+i+NUM_DIRECTIONS)%NUM_DIRECTIONS;
        si=x+dx[dd];
        sj=y+dy[dd];
		if((si>=0)&&(si<w)&&(sj>=0)&&(sj<h)) {
        	if(ras_edgemap(sj,si) == 2) return;
		}
    } /* not blocked */

    cl = ras_edgemap(y,x);
    for(i=le;i<=ue;i++)
	{
        dd=(d+i+NUM_DIRECTIONS)%NUM_DIRECTIONS;
        si=x+dx[dd];
        sj=y+dy[dd];
		if((si>=0)&&(si<w)&&(sj>=0)&&(sj<h))
		{
     	   	if(ras_edgemap(sj,si) == 1)
			{
       	     	ras_edgemap(sj,si) = cl;
       	     	IlePropagate(ras_edgemap,si,sj,dd);
       	 	}
		}
    }
}

void MCMC_EdgeLinking::DetectedEdgesToImage(void)
{
	int		k,m,i,j;
	MCMC_OpAnEdge *pedge;

	for (k=0; k<NumEdges(); k++)
	{
		pedge = &detected_edges.Edge(k);
		for (m=0; m<pedge->NumPoints(); m++)
		{
			j = pedge->j(m);
			i = pedge->i(m);

			ras_edges(j, i) = k+3;
		}
	}
}

void MCMC_EdgeLinking::ConnectEdges(void)
//iptr;	edge map
//eptr;	edge image
{
	// This routine assigns the junction pixels. It 				
	// assigns junction points to edges that touch them            	
	// it also forms the junction data base, and sets junctions    	
	// that have been processed to 0. Finally, it checks to see if 	
	// any junctions are left, and prints an error if there are    

	// modifications:	 do not assign junction pixels to lines!

	int w,h,i,j;
	int num_junctions;

	int ii,jj;
	int cl,assigned,num,numlist;
	w = cols();
	h = rows();

	num_junctions=0;
	for(i=1;i<(w-1);i++) {
		for(j=1;j<(h-1);j++) {
			if(ras_edges(j,i)==2) num_junctions++;
		}
	}

	mcmcVector<int> vc_tlist;

	vc_tlist.SetRowVector(NumEdges());

	num_junctions=0;
	for(i=1;i<(w-1);i++) {
		for(j=1;j<(h-1);j++) {
			if(ras_edges(j,i)==2) { // junction
				vc_tlist.InitValue(0);
				numlist=0;

//ztu			ptr=jptr->junction[num_junctions];
// Modification by Brendan McCane Nov 4, 1993
//ztu           ptr->x_pos = i;
//              ptr->y_pos = j;
// End mod
				num_junctions++;
				assigned=0;
				for(ii= -1;ii<=1;ii++) {
					for(jj= -1;jj<=1;jj++) {
						if((ii==0)&&(jj==0)) continue;
						cl=ras_edges(j+jj,i+ii);
						if(cl>=3) { // found a valid edge neighbour
							num=cl-3;
							if(vc_tlist(num)==0) { // not assigned to a given edge yet
								detected_edges.Edge(num).InsertPoint(j, i);
								assigned=1;
								// assign to junction
								vc_tlist(num)=1;
								numlist++;
							}
						}
					}
				}
				if(assigned)
					ras_edges(j,i)=0;
			}
		}
	}

	detected_edges.CopyFromPixelList();
}

int MCMC_EdgeLinking::NumNeighbors(Raster<int> &ras_image, int x, int y, int neighbourhood)
{
	int ii, jj;
	int num_neighbours=0;
	int width, height;
	int x_start, y_start, x_fin, y_fin;

	width = ras_image.cols();
	height = ras_image.rows();
	if (x-neighbourhood>=0) x_start = x-neighbourhood;
	else x_start=0;
	if (x+neighbourhood<width) x_fin = x+neighbourhood;
	else x_fin = width-1;
	if (y-neighbourhood>=0) y_start = y-neighbourhood;
	else y_start=0;
	if (y+neighbourhood<height) y_fin = y+neighbourhood;
	else y_fin = height-1;

	for (ii=x_start; ii<=x_fin; ii++)
		for (jj=y_start; jj<=y_fin; jj++)
			if (ras_image(jj, ii))
				num_neighbours++;

	return(num_neighbours);
}

void MCMC_EdgeLinking::BrenDda(MCMC_OpAnEdge& edge_seg, int x1, int y1, int x2, int y2)
{
	edge_seg.Empty();

	int dx, dy, steps, i;
	float x_inc, y_inc, x, y;

	dx = x2 - x1;
	dy = y2 - y1;
	if (abs(dx)>abs(dy)) steps = abs(dx);
	else steps = abs(dy);

	if (steps==0)
		return;

	x_inc = (float) dx/steps;
	y_inc = (float) dy/steps;
	x = (float)x1;
	y = (float)y1;

	int m,n;

	for (i=0; i<=steps; i++)
	{
		n  = (int)(x+0.5);
		m = (int)(y+0.5);
		
		edge_seg.InsertPoint(m, n);

		x += x_inc;
		y += y_inc;
	}

	edge_seg.CopyFromPixelList();
}

void MCMC_EdgeLinking::ExtendEdgeNew(MCMC_OpAnEdge &edge_new, MCMC_OpAnEdge &edge_old, int end, int smooth, int ext)
{
	int delta_x, delta_y, newx, newy, i;
	
	if (edge_old.NumPoints()<=smooth)
		smooth = edge_old.NumPoints()-1;

	if (smooth<1)
		return;

	/* calculate the gradient of the edge */
	int m1,n1,m2,n2;

	i = 0;
	m1 = edge_old.j(end);
	n1 = edge_old.i(end);

	if (end==0)
	{
		m2 = edge_old.j(end+smooth);
		n2 = edge_old.i(end+smooth);
	}
	else
	{
		m2 = edge_old.j(end-smooth);
		n2 = edge_old.i(end-smooth);
	}
	delta_x = n1 - n2;
	delta_y = m1 - m2;

	/*
	 * set the new points to be extrapolated proportional to the amount 
	 * of extension.
	 */
	if (ext>edge_old.NumPoints()) ext = edge_old.NumPoints();
	newx = n1+delta_x*ext/smooth;
	newy = m1+delta_y*ext/smooth;

	/* form the new edge */
	BrenDda(edge_new, n1, m1, newx, newy);
	edge_new.DeleteTheFirstPixel();
}

#define bren_in_interval(a, start, end) 						\
					(((((a)>=(start))&&((a)<=(end)))||		\
				 	(((a)>=(end))&&((a)<=(start))))?1:0)


float euclid_dist(int x1, int y1, int x2, int y2)
{
	float ans;

	ans = (float)sqrt((double)(x1-x2)*(double)(x1-x2)+
					   (double)(y1-y2)*(double)(y1-y2));
	return(ans);
}

#define END_EDGE 20000
#define START_EDGE 10000
void MCMC_EdgeLinking::handle_runedge(int run_edge_no, int cx, int cy,
									  MCMC_OpEdges *start_extend, MCMC_OpEdges *end_extend, MCMC_OpEdges *start_use,
									  MCMC_OpEdges *end_use)
{
	MCMC_OpAnEdge *run_edge, new_edge;
	int		i, index=0;
	float	dist, min_dist=1e10;
	int		width, height;

	width = cols();
	height = rows();
	if (run_edge_no>=END_EDGE)
	{
		run_edge = &end_extend->Edge(run_edge_no-END_EDGE);
		/* we ran into an extended edge */
		if (run_edge->NumPoints()>0)
		{
			/*
			 * cut the extended edge back to the junction pt
			 * and add to the end database.
			 */
			for (i=0; i<run_edge->NumPoints(); i++)
			{
				if ((dist=euclid_dist(cx, cy, run_edge->x(i),
								run_edge->y(i)))<2.9)
				{
					if (dist<min_dist)
					{
						min_dist = dist;
						index = i;
					}
				}
			}
			new_edge.CopyFrom(*run_edge, index+1);
			
			if (min_dist>1.5)
			{
				int addx, addy;
				addx=(int)((float)(cx+run_edge->x(new_edge.NumPoints()-1))/2.0);
				addy=(int)((float)(cy+run_edge->y(new_edge.NumPoints()-1))/2.0);
				new_edge.AddPoint(addy, addx);
			}

			end_use->Edge(run_edge_no-END_EDGE) = new_edge;
			for (i=0; i<run_edge->NumPoints(); i++)
			{
				if ((bren_in_interval(run_edge->x(i), 0, width-1))&&
					(bren_in_interval(run_edge->y(i), 0, height-1))&&
					(ras_edges(run_edge->y(i), run_edge->x(i)) == run_edge_no))
					ras_edges(run_edge->y(i), run_edge->x(i)) = 0;
			}
			for (i=0; i<new_edge.NumPoints(); i++)
			{
				if ((bren_in_interval(new_edge.x(i), 0, width-1))&&
					(bren_in_interval(new_edge.y(i), 0, height-1)))
					ras_edges(new_edge.y(i), new_edge.x(i)) =	run_edge_no;
			}
			end_extend->Edge(run_edge_no-END_EDGE).Empty();
		}
	}
	else if (run_edge_no>=START_EDGE)
	{
		run_edge = &start_extend->Edge(run_edge_no-START_EDGE);
		/* we ran into an extended edge */
		if (run_edge->NumPoints()>0)
		{
			/*
			 * cut the extended edge back to the junction pt
			 * and add to the start database.
			 */
			for (i=0; i<run_edge->NumPoints(); i++)
			{
				if ((dist=euclid_dist(cx, cy, run_edge->x(i),
								run_edge->y(i))<2.9))
				{
					if (dist<min_dist) min_dist = dist;
					index = i;
				}
			}
			new_edge.CopyFrom(*run_edge, index+1);

			if (min_dist>1.5)
			{
				int addx, addy;
				addx=(int)((float)(cx+run_edge->x(new_edge.NumPoints()-1))/2.0);
				addy=(int)((float)(cy+run_edge->y(new_edge.NumPoints()-1))/2.0);
				new_edge.AddPoint(addy, addx);
			}
			start_use->Edge(run_edge_no-START_EDGE) = new_edge;
			for (i=0; i<run_edge->NumPoints(); i++)
			{
				if ((bren_in_interval(run_edge->x(i), 0, width-1))&&
					(bren_in_interval(run_edge->y(i), 0, height-1))&&
					(ras_edges(run_edge->y(i), run_edge->x(i)) == run_edge_no))
					ras_edges(run_edge->y(i), run_edge->x(i)) = 0;
			}
			for (i=0; i<new_edge.NumPoints(); i++)
			{
				if ((bren_in_interval(new_edge.x(i), 0, width-1))&&
					(bren_in_interval(new_edge.y(i), 0, height-1)))
					ras_edges(run_edge->y(i), run_edge->x(i))=	run_edge_no;
			}
			start_extend->Edge(run_edge_no-START_EDGE).Empty();
		}
	} 
}


void MCMC_EdgeLinking::ExtendParallelEdges(MCMC_OpEdges &local_edges, int smooth, int max_ext)
{
	int time, edge_no;
	MCMC_OpEdges start_extend, end_extend, start_use, end_use;
	int width, height;
	int run_edge_no,i,j;

	local_edges = detected_edges;

	width = cols();
	height = rows();

	start_use.SetNumEdges(detected_edges.NumEdges());
	end_use.SetNumEdges(detected_edges.NumEdges());
	start_extend.SetNumEdges(detected_edges.NumEdges());
	end_extend.SetNumEdges(detected_edges.NumEdges());

	MCMC_OpAnEdge *edge,*end_edge;

	for (edge_no=0; edge_no<detected_edges.NumEdges(); edge_no++)
	{
		edge = &detected_edges.Edge(edge_no);
		assert(edge->GetFirst(j, i));

		if (NumNeighbors(ras_edges, i, j, 1)==2)
			ExtendEdgeNew(start_extend.Edge(edge_no), *edge, 0, smooth, max_ext);;

		assert(edge->GetLast(j, i));
		if (NumNeighbors(ras_edges, i, j, 1)==2)
			ExtendEdgeNew(end_extend.Edge(edge_no), *edge, edge->NumPoints()-1, smooth, max_ext);
	}
	
	for (time=0; time<max_ext; time++)
	{
		for (edge_no=0; edge_no<detected_edges.NumEdges(); edge_no++)
		{
			int cx, cy, i, j;
			int foundx, foundy, newx, newy;
			if (end_extend.Edge(edge_no).NumPoints() > 0)
			{
				end_edge = &end_extend.Edge(edge_no);

				if (end_edge->NumPoints()<=time) continue;
				
				if ((bren_in_interval(end_edge->x(time), 0, width-1))&&
					(bren_in_interval(end_edge->y(time), 0, height-1)))
					ras_edges(end_edge->y(time), end_edge->x(time)) =
						edge_no+END_EDGE;
				/*
				 * if the number of neighbours > 2 stop add the current
				 * edge to the use_edge database and also remove the edge it
				 * ran into if appropriate.
				 */
				if (NumNeighbors(ras_edges, end_edge->x(time), end_edge->y(time), 2)>3)
				{
					cx = end_edge->x(time);
					cy = end_edge->y(time);
					/* create the edge to add */
					end_use.Edge(edge_no).CopyFrom(*end_edge, time+1);

					/* free the extended edge */
					end_extend.Edge(edge_no).Empty();

					/* track down the edge we ran into */
					run_edge_no = -1;
					for (i=-2; i<=2; i++)
					{
						for (j=-2; j<=2; j++)
						{
							if ((bren_in_interval(cx+i, 0, width-1))&&
								(bren_in_interval(cy+j, 0, height-1))&&
								ras_edges(cy+j, cx+i) &&
								(ras_edges(cy+j, cx+i)-END_EDGE != edge_no) &&
								(ras_edges(cy+j, cx+i)-3!=edge_no))
							{
								run_edge_no = ras_edges(cy+j, cx+i);
								foundx = cx+i;
								foundy = cy+j;
							}
						}
					}
					if (run_edge_no>=0)
					{
						newx = (int)((cx+foundx)/2.0);
						newy = (int)((cy+foundy)/2.0);
						end_use.Edge(edge_no).AddPoint(newy, newx);

						handle_runedge(run_edge_no, newx, newy, &start_extend,
									   &end_extend, &start_use, &end_use);
					}
				}
				else if ((end_edge->x(time)==0)||(end_edge->x(time)==width)||
						 (end_edge->y(time)==0)||(end_edge->y(time)==height))
				{
					cx = end_edge->x(time);
					cy = end_edge->y(time);

					end_use.Edge(edge_no).CopyFrom(*end_edge, time+1);

					end_extend.Edge(edge_no).Empty();
				}
			}
			if (start_extend.Edge(edge_no).NumPoints() > 0)
			{
				MCMC_OpAnEdge *start_edge=&start_extend.Edge(edge_no);

				if (start_edge->NumPoints()<=time) continue;
				if ((bren_in_interval(start_edge->x(time), 0, width-1))&&
					(bren_in_interval(start_edge->y(time), 0, height-1)))
					ras_edges(start_edge->y(time), start_edge->x(time))
						= edge_no+START_EDGE;
				/*
				 * if the number of neighbours > 2 stop add the current
				 * edge to the use_edge database and also remove the edge it
				 * ran into if appropriate.
				 */
				if (NumNeighbors(ras_edges, start_edge->x(time), start_edge->y(time), 2)>3)
				{
					int i;

					cx = start_edge->x(time);
					cy = start_edge->y(time);
					/* create the edge to add */
					start_use.Edge(edge_no).CopyFrom(*start_edge, time+1);


					/* free the extended edge */
					start_extend.Edge(edge_no).Empty();

					/* track down the edge we ran into */
					run_edge_no = -1;
					for (i=-2; i<=2; i++)
					{
						for (j=-2; j<=2; j++)
						{
							if ((bren_in_interval(cx+i, 0, width-1))&&
								(bren_in_interval(cy+j, 0, height-1))&&
								(ras_edges(cy+j, cx+i))&&
								(ras_edges(cy+j, cx+i)-START_EDGE!=edge_no)&&
								(ras_edges(cy+j, cx+i)-3!=edge_no))
							{
								run_edge_no=ras_edges(cy+j, cx+i);
								foundx = cx+i;
								foundy = cy+j;
							}
						}
                                                                                
					}
					if (run_edge_no>=0)
					{
						newx = (int)((cx+foundx)/2.0);
						newy = (int)((cy+foundy)/2.0);
						start_use.Edge(edge_no).AddPoint(newy, newx);
						handle_runedge(run_edge_no, newx, newy, &start_extend, 
								   &end_extend, &start_use, &end_use);
					}
				}
				else if ((start_edge->x(time)==0)||
						 (start_edge->x(time)==width)||
						 (start_edge->y(time)==0)||
						 (start_edge->y(time)==height))
				{
					cx = start_edge->x(time);
					cy = start_edge->y(time);

					start_use.Edge(edge_no).CopyFrom(*start_edge, time+1);

					/* free the extended edge */
					start_extend.Edge(edge_no).Empty();
				}
			}
		}
	}

	for (edge_no=0; edge_no<local_edges.NumEdges(); edge_no++)
	{
		char change;
		if (end_use.Edge(edge_no).NumPoints() > 0)
		{
			MCMC_OpAnEdge *edge=&end_use.Edge(edge_no);
			for (i=0; i<edge->NumPoints(); i++)
				local_edges.Edge(edge_no).InsertPoint(edge->y(i), edge->x(i));
										   
			change=true;
		}
		if (start_use.Edge(edge_no).NumPoints() > 0)
		{
			MCMC_OpAnEdge *edge=&start_use.Edge(edge_no);
			for (i=0; i<edge->NumPoints(); i++)
				local_edges.Edge(edge_no).InsertPoint(edge->y(i), edge->x(i));

			change=true;
		}
		local_edges.Edge(edge_no).CopyFromPixelList();
	}
}

// Label number on edge pixels that form connected line
void MCMC_EdgeLinking::LabelNumberOnEdgePixels(Raster<int> &ras_edgemap)
{
	int	w,h,counte,countj,i,j,ii,jj;
	int numsegments,count;
	int si,sj,done,d,pass,found;

	w = ras_edgemap.cols();
	h = ras_edgemap.rows();

	/* do a little bit of error checking */

	for(i=1; i<(w-1); i++)
	{
		for(j=1; j<(h-1); j++)
		{
			if(ras_edgemap(j,i) == 1)
			{	/* edge pixel */
				counte=0;
				countj=0;
				for(ii=-1; ii<=1; ii++)
				{
					for(jj=-1; jj<=1; jj++)
					{
						if((ii==0) && (jj==0)) continue;
						if(ras_edgemap(j+jj,i+ii) == 1) counte++;
						if(ras_edgemap(j+jj,i+ii) == 2) countj++;
					}
				}
			}
		}
	}

	/* junctions are 2, edges are 1 */

	numsegments=0;
	done=0;
	while(!done)
	{
		done=1;
		found=0;
		for(i=1;(i<(w-1))&&(!found);i++)
		{
			for(j=1;(j<(h-1))&&(!found);j++)
			{
				if(ras_edgemap(j,i) == 1)
				{
					/* have to check neighbourhood for junctions */
					found=0;
					count=0;
					pass=1;
					si=i;
					sj=j;
					for(ii=-1; ii<=1; ii++)
					{
						for(jj=-1; jj<=1; jj++)
						{
							if((ii==0) && (jj==0)) continue;
							if(ras_edgemap(j+jj,i+ii) == 2) pass=0; /* dont use */
							if(ras_edgemap(j+jj,i+ii) == 1) count++;
						}
					}
					if((pass)&&(count==2)) found=1;
				}
			}
		}
		if(found)
		{
			ras_edgemap(sj,si) = numsegments+3;
			for(d=0; d<NUM_DIRECTIONS; d++)
				IlePropagate(ras_edgemap,si,sj,d);
			done=0;
			numsegments++;
		}
	}

	/* for closed, no-junctions regions, start anywhere */
    done=0;
    while(!done)
	{
        done=1;
        found=0;
        for(i=1;(i<(w-1))&&(!found);i++)
		{
            for(j=1;(j<(h-1))&&(!found);j++)
			{
				if (ras_edgemap(j,i) == 1)
				{
                    /* have to check neighbourhood for junctions */
                    found=0;
                    count=0;
                    pass=1;
                    si=i;
                    sj=j;
                    for(ii=-1; ii<=1; ii++)
					{
                        for(jj=-1; jj<=1; jj++)
						{
                            if((ii==0) && (jj==0)) continue;
                            if(ras_edgemap(j+jj,i+ii) == 2) pass=0; /* use */
                            if(ras_edgemap(j+jj,i+ii) == 1) count++;
                        }
                    }
                    if((pass==1)&&((count==1)||(count==2))) found=1;
                }
            }
        }
        if(found)
		{	/* texture line or int edge */
			ras_edgemap(sj,si) = numsegments+3;
            for(d=0; d<NUM_DIRECTIONS; d++)
                IlePropagate(ras_edgemap,si,sj,d);
            done=0;
            numsegments++;
        }
    }

	/* check that no 1 pixels remain */
	for(i=1;i<(w-1);i++)
	{
		for(j=1;j<(h-1);j++)
		{
			if(ras_edgemap(j,i) == 1)
			{
/*
 * Modification by Brendan McCane to include isolated pixels as an edge
 * and edges with only two pixels in them.
 */
                ras_edgemap(j,i) = numsegments+3;
                for (ii=-1; ii<=1; ii++)
                for (jj=-1; jj<=1; jj++)
                {
                    if ((ii==0) && (jj==0)) continue;
                    if (ras_edgemap(j+jj,i+ii) == 1)
                        ras_edgemap(j+jj,i+ii) = ras_edgemap(j,i);
                }
                numsegments++;
			}
		}
	}
}

#define NUM_DIRECTIONS 8
